1.2.1 사용한 패키지

자료손질시 사용하였던 패키지 목록입니다.

# 패키지 로드
library(readxl)  #엑셀형식
library(tidyverse)  #이제는 일단깔고
library(plyr)  # 데이터프레임 다룰 때 사용

1.2.2 자료 불러오기

반복-조건문을 사용하여 12개월치 자료를 순차적으로 불러온 뒤, 한 데이터프레임으로 만들었습니다.

## data de air
fi <- paste0("./data/2019AQdata/", list.files(path = "./data/2019AQdata/", pattern = "확정자료엑셀다운*"))

for (i in 1:12) {
    assign(paste0("A", i, "raw"), read_xls(fi[i]))
    if (i == 1) {
        Araw <- plyr::rbind.fill(A1raw)
        rm(A1raw)  # 부산물 삭제
    } else {
        Araw <- plyr::rbind.fill(Araw, eval(as.symbol(paste0("A", i, "raw"))))
        rm(list = as.character(paste0("A", i, "raw")))  #부산물 삭제
    }
    
}

Araw

           날짜      시도 측정소명 측정소코드 아황산가스 일산화탄소 오존
1 2019-01-01 01 서울 중구     중구     111121       .003         .7 .003
2 2019-01-01 02 서울 중구     중구     111121       .003         .8 .002
3 2019-01-01 03 서울 중구     중구     111121       .003         .9 .002
4 2019-01-01 04 서울 중구     중구     111121       .003         .8 .002
5 2019-01-01 05 서울 중구     중구     111121       .003         .8 .002
6 2019-01-01 06 서울 중구     중구     111121       .003         .8 .002
7 2019-01-01 07 서울 중구     중구     111121       .003         .7 .006
  이산화질소 PM10 PM2.5
1       .054   39    24
2       .056   38    27
3       .057   42    28
4       .054   42    31
5       .048   49    33
6       .048   47    32
7        .04   48    30
 [ reached 'max' / getOption("max.print") -- omitted 218993 rows ]

1.2.3 자료의 일 평균 구하기

시간대별로 수치가 기록되어있는 자료를 일평균을 내주는 작업입니다.

먼저 측정소 순으로 정렬한 뒤, YYYY-MM-DD HH로 이루어진 날짜 형식을 날짜와 시간으로 분리하여 일평균을 구할 수 있도록 날짜 열을 만들었습니다. 1

Araw$PM10 <- as.numeric(Araw$PM10)

order_Air <- Araw[c(order(Araw$측정소코드)), ]


Araw <- separate(data = Araw, col = 날짜, into = c("날짜", "시"), sep = " ")

다음은 일 평균을 구하는 과정입니다.

davg <- vector(mode = "numeric", length(Araw$측정소코드)/24)

for (j in 1:as.numeric(length(Araw$측정소코드)/24)) {
    tmp <- vector(mode = "numeric", 24)
    for (k in 1:24) {
        if (is.null(Araw$PM10[(j - 1) * 24 + k]) == TRUE) {
            davg[j] <- NULL
            break
        } else {
            tmp[k] <- Araw$PM10[(j - 1) * 24 + k]
        }
    }
    davg[j] <- mean(tmp)
}

rm(tmp)

dmean_air <- (Araw[c(seq(from = 1, to = length(Araw$PM10), by = 24)), ] %>% cbind(davg))[, 
    c(1, 3:5, 12)]

dmean_air

          날짜      시도 측정소명 측정소코드      davg
1   2019-01-01 서울 중구     중구     111121  35.95833
25  2019-01-02 서울 중구     중구     111121  33.04167
49  2019-01-03 서울 중구     중구     111121  35.25000
73  2019-01-04 서울 중구     중구     111121  54.70833
97  2019-01-05 서울 중구     중구     111121  60.04167
121 2019-01-06 서울 중구     중구     111121  42.41667
145 2019-01-07 서울 중구     중구     111121  50.75000
169 2019-01-08 서울 중구     중구     111121  37.66667
193 2019-01-09 서울 중구     중구     111121  46.62500
217 2019-01-10 서울 중구     중구     111121  51.25000
241 2019-01-11 서울 중구     중구     111121  70.95833
265 2019-01-12 서울 중구     중구     111121  89.04167
289 2019-01-13 서울 중구     중구     111121 106.66667
313 2019-01-14 서울 중구     중구     111121 155.20833
337 2019-01-15 서울 중구     중구     111121 122.45833
 [ reached 'max' / getOption("max.print") -- omitted 9110 rows ]

위처럼 반복문을 사용하여, 24개의 열씩 합하였습니다.

1.2.4 시간대 추출

해당 분석을 위하여 19년도 미세먼지를 검색하던 도중 눈에 띄는 뉴스를 하나 보게 되었습니다.

위 뉴스를 접하고 1월중 가장 미세먼지가 심한 날을 알아보았습니다.

max <- 0  #init value
date <- "2019-01-01"

for (j in 1:31) {
    # 31일까지
    if (isTRUE(dmean_air$davg[j] > max)) {
        max <- dmean_air$davg[j]
        date <- dmean_air$날짜[j]
    } else {
        max <- max
        date <- date
    }
}

max

[1] 155.2083

date

[1] "2019-01-14"

신기하게도 뉴스가 방영된 날짜와 똑같음을 알 수 있었습니다.

해당 사례로 분석을 진행하기 위해 19년 1월 14일의 관측소 데이터만 추출하였습니다.

maxair <- dmean_air[0, ]

빈 데이터프레임을 만들고,

for (j in 1:length(dmean_air$davg)) {
    if (dmean_air$날짜[j] == date) {
        maxair[j, ] <- dmean_air[j, ]
    }
}

maxpm <- filter(maxair, !is.na(davg))

반복문을 사용하여 추출한 뒤, 결측치를 제거하는 방식으로 빈 열들을 제거하였습니다.

1.2.5 기존 자료와 주소 조인

해당 관측소의 위치를 알기 위하여 에어코리아에서 구득한 관측소 위치자료를 기존 데이터프레임에 조인하였습니다.

lctn <- read_xls("./data/2019AQdata/station_list.xls", range = "B4:C29")

maxpm_refined <- inner_join(maxpm, lctn, by = "측정소명")

1.3.1 사용한 패키지

네이버 클라우드에서 제공하는 지오코딩 API 서비스를 사용하기 위한 API 관련 패키지입니다. 기존 구글로 지오코딩하는 방식은 지도의 최신화가 안되는 등의 이유로 국내 포털보다 정확도가 떨어져서 네이버 지오코딩 서비스를 사용하였습니다. 네이버 클라우드

library(httr)
library(XML)

1.3.2 지오코딩

우선, 발급받은 API ID와 비밀번호를 입력하였습니다. (API키는 개인 키라 가렸습니다.)

# Client ID 설정
clientId <- "API ID"

# Client Secret 설정
clientSecret <- "API PSSWD"

API 호출을 통해 지오코딩을 하는 과정입니다. 기존 API 호출로는 하나의 주소밖에 처리하지 못하므로, 이 또한 반복문을 사용하여 한번에 처리될 수 있게 하였습니다.

x <- vector(mode = "double", length(maxpm_refined$davg))
y <- vector(mode = "double", length(maxpm_refined$davg))

for (l in 1:length(maxpm_refined$davg)) {
    
    apiResult <- httr::GET(url = "https://naveropenapi.apigw.ntruss.com/map-geocode/v2/geocode", 
        httr::add_headers(`X-NCP-APIGW-API-KEY-ID` = clientId, `X-NCP-APIGW-API-KEY` = clientSecret), 
        query = list(query = maxpm_refined$`측정소 주소`[l]  #주소
)  #리버스지오코딩시 쿼리 추가
)
    if (apiResult$status_code == "200") {
        result <- rawToChar(apiResult$content)
        Encoding(result) <- "UTF-8"
        list <- xmlToList(result)
        
        u <- as.numeric(list$addresses$x)
        v <- as.numeric(list$addresses$y)
        
        x[l] <- u
        y[l] <- v
        
        
    } else if (apiResult$status_code == "400") {
        print(paste("Bad Request Exception in", maxpm_refined$측정소명[l]))
        
    } else if (apiResult$status_code == "500") {
        print(paste("Unexpected Error in", maxpm_refined$측정소명[l]))
    } else {
        print("몰라요")
    }
    
}

maxpm_refined$x <- x
maxpm_refined$y <- y

1.4.1 사용한 패키지

크리깅 보간 및 도식을 위해 사용한 패키지입니다.

library(sf)
library(tmap)
library(rgdal)
library(raster)
library(leafem)

1.4.2 서울시 행정경계

eTL에서 구득한 서울시 행정경계 파일을 마스킹을 통한 추출 및 CRS획득을 위해 불러왔습니다.

seoulpoly <- readOGR(dsn = "./data/shp", layer = "tl_emd._seoul_4326", encoding = "UTF-8", 
    stringsAsFactors = FALSE)

OGR data source with driver: ESRI Shapefile 
Source: "C:\Users\dydgn\OneDrive\바탕 화면\Analytical_Methods_for_Spatial_Information_Assignment\ASS-11\data\shp", layer: "tl_emd._seoul_4326"
with 467 features
It has 3 fields

crs <- crs(seoulpoly)

1.4.3 교차검증을 위한 Samples

3개의 관측소 데이터를 무작위로 추출한 후, 분석하려는 데이터프레임에서 해당 자료를 제거하였습니다. 추출한 데이터는 추후 교차검증시 사용하였습니다.

PMdataCV <- sample_n(maxpm_refined, 3)
PMdataCVpts <- SpatialPointsDataFrame(PMdataCV[, 7:8], PMdataCV, proj4string = crs)


desmpld <- maxpm_refined[!(maxpm_refined$davg == PMdataCV$davg[1] | maxpm_refined$davg == 
    PMdataCV$davg[2] | maxpm_refined$davg == PMdataCV$davg[3]), ]

1.4.4 공간자료로의 변환

지오코딩하여 좌표를 획득한 자료를 공간자료로 변환하였습니다.

PMdatapts <- SpatialPointsDataFrame(desmpld[, 7:8], desmpld, proj4string = crs)

tmap_mode("view")
tm1 <- tm_shape(PMdatapts) + tm_dots("davg", title = "일평균 PM10 <br/>(㎍/m³)", 
    id = "측정소명", legend.show = TRUE)
lf <- tmap_leaflet(tm1)
addMouseCoordinates(lf)